arXiv:1506.03848vl [math.NA] llJun2015 


On evaluation of the Heun functions 


Oleg V. Motygin 

Institute of Problems in Mechanical Engineering, Russian Academy 
of Sciences, V.O., Bol’shoy pr., 61, 199178 St. Petersburg, Russia 
email: o.v.motygin@gmail.com 


Abstract 

In the paper we deal with the Heun functions — solutions of the Heun equation, which is the 
most general Fuclrsian equation of second order with four regular singular points. Despite the 
increasing interest to the equation and numerous applications of the functions in a wide variety 
of physical problems, it is only Maple amidst known software packages which is able to evaluate 
the Heun functions numerically. But the Maple routine is known to be imperfect: even at regular 
points it may return infinities or end up with no result. Improving the situation is difficult because 
the code is not publicly available. The purpose of the work is to suggest and develop alternative 
algorithms for numerical evaluation of the Heun functions. A procedure based on power series 
expansions and analytic continuation is suggested which allows us to avoid numerical integration 
of the differential equation and to achieve reasonable efficiency and accuracy. Results of numerical 
tests are given. 


1 Introduction 

In the present paper we deal with the Heun functions which are solutions of the equation introduced 
by Karl Heun in 1889 [3] as a generalization of the hypergeometric equation. Heun’s equation is 
the most general Fuclrsian equation of second order with four regular singular points; we refer to 
[12, 14, 9, 15] for a comprehensive mathematical treatment of the topic. At the same time, it is 
only in recent years when the equation has become popular in the physics literature. Now the Heun 
equation appears in many fields of modern physics and it is used to describe a wide variety of physical 
phenomena — a comprehensive literature on physical applications can be found in [4]. As a good 
source of information on the current development of the held one should also mention “The Heun 
project”: http://theheunproject.org/. 

Despite the increasing interest to the equation, at the moment the only, to author’s knowledge, 
software package able to evaluate the Heun functions numerically is Maple. However, the code is far 
from being perfect, it is not difficult to encounter problems, when for ordinary parameters, at regular 
points the routine returns infinities or spends tens seconds with no result. The quality of the code is a 
serious obstacle for (potentially very numerous) applications of the functions. Improving the situation 
is virtually impossible because the code is not publicly available. 

The purpose of the present work is to develop alternative algorithms for numerical evaluation of the 
Heun functions. We suggest a procedure based on power series expansions and analytic continuation 
which allows us to avoid numerical integration of the differential equation and to achieve reasonable 
efficiency and accuracy. Program code is presented in [11]. Results of numerical tests and comparison 
with a case when Heun’s function reduces to a simple algebraic function are given. 

The algorithm is applicable for computation of the multi-valued Heun functions. In the present 
paper we also define their single-valued counterparts by fixation of branch cuts. (Notably, in most 
studies of the Heun functions the subject is not discussed; the author has also been unable to find 
information on branch cuts in Maple’s documentation on the Heun function.) For the single-valued 
functions an improvement of the algorithm for points close to the singular ones is described. 

It should be noted that the developed algorithms are not intended to be universal — more or less 
ordinary parameters are assumed. Surely, numerical problems are expected and special treatment is 
needed for the cases of merging singular points (see [7]) or large accessory parameter. 
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2 Statement and basic notations 


We start by writing Heun’s equation in the standard form (see [2]) 
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The parameter g £ C is usually referred to as an accessory or auxiliary parameter and a, (5, 7 , S, 
e (also belonging to C) are exponent-related parameters connected via the Fuchsian relation 


a P P 1 — 7 + <5 + e. 


The equation has four regular singular points located at z = 0. 1, a, 00 . It will be assumed below that 
a G C, a / {0,1, 00 }. In the notation of the solutions sometimes we will omit some of the parameters 
so that H(z) = H{a,... \ z) = H(a, q, a, (3, 7 , 5; z). 

The Frobenius method can be used to derive local power-series solutions to (1). There are 8 local 
solutions of equation (1) (two per a singular point). In § 3 we will present the two local power-series 
solutions in a neighbourhood of the point z = 0. One of them is analytic in a vicinity of zero and if 7 
is not a nonpositive integer, this solution, normalized to unity, is called the local Heun function (see 
[12]). It is usually denoted by HI (a, q, a, f3, 7 , <5; z). For the second Frobenius local solution we will use 
the notation Hs(a, q, a, /3, 7 , <5; z). 

When 7 is a nonpositive integer, one solution of (1) is analytic but equal to zero at z = 0, whereas 
the second solution can be normalized to unity but generally is not analytic. It can be an arguable 
point, but we will denote by HI (a, q, a, /?, 7 , 5] z) the normalized solution and by Hs(a, q, a, /3, 7 , <5; z) 
the analytic one. 

It is known that generally HI (a,. ..; z) is a multi-valued function with branch points at z = 1, a 
and 00 and, so, to define a single-valued function we should choose branch cuts. In the present work 
we fix the branch cuts &ioo = (l,+oo) and SS a oo = (a, e iarg ^oo) connecting the points 1 and a to 
00 , respectively (see fig. 1). The second function Hs(a,...; z ) has the same branch cuts &ioc and 
SSaoo- Besides, the function Hs(z) — generally, and the function Hl(z) — for 7 e {0}UZ' (Z!~ means 
the set of negative integers), have a branch point at z = 0. So, we will also define the branch cut 
&O 00 = (-oo,0). 

It should be noted that the choice of branch cuts in Maple is unclear, but, definitely, differs from 
that in the present paper. With the present fixation the function definition domain is star-like with 
respect to zero which is natural as we will define functions Hl(z), Hs(z) by analytic continuation from 
a vicinity of z = 0 (see § 5). Besides, there is good agreement between branch cuts of the single-valued 
functions HI, Hs and their representations near singular points (see § 6 ). 
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3 Power series expansion at the point z = 0 

Power series expansion of the local Heun function HI, such that Hl(a, q, a, (5, 7 , h; 0) = 1, is well-known 
since [3] for 7 7 {0} U Z _ . We have 

OO 

Hl(a, q , a, (5, 7 , 6;z) = Y b n z n , (2) 

71=0 

where the coefficients b n are subject to the following three-term recurrence: 

Pnbn = Qnbn— 1 “I - Rnbn— 2 - (3) 


Here 

P n = an (7 - 1 + n), 

Qn = q + (n - 1 ) [(a + 1)(7 + n — 2 ) + e + a<5], 
R n = — (n — 2 + a){n — 2 + /?) 

and the initial conditions are 

b~ i = 0, bo = 1. 


(4) 


The Heun function Hl(z ) is analytic in the circle |z| < i? 0 i where 7?o is the distance from zero to the 
nearest singular point, i?o = min{l, |a|}. Then the famous Cauchy’s theorem on the expansion of an 
analytic function into a power series (see e.g. Theorem 16.7 in [10, Part 1]) guarantees that the series 
( 2 ) converges to HI inside the circle \z\ < Rq- At this, of course, there might be issues with stability of 
the recurrence process. For the stability it can be useful to write (2) in the form Hl(z ) = o b n (z), 
where b n (z) = b n z n and the recurrence takes the form: P n b n = zQ n b n -\ + z 2 R n b n - 2 . 

In the case the local Frobenius solution corresponding to the smaller exponent (0 or 1 — 7 ) 

may contain a logarithmic factor (see e.g. [5]). So, for 7 £ {0} UZ” we are looking for the solution of 
(1), satisfying Hl( 0) = 1, in the following form: 


HI (a, q, a, (5, 7 , <5; z) 


OO OO 

^2 c nZ H + log(z) ^ 

n=0, 717^71* 71=71* 


(5) 


where 71 * = 1 — 7 . Note that a solution with the sought property Hl( 0) = 1 could be found for any 
c nt . We fix c n< , = 0 for definiteness. 

Substituting (5) into (1) we find 


log (z)Jfl 


SnZ 


+ J?I Y CnZ n ) +J?I Y S n Zn j = 0 ’ 

\n=0,n^n* / \ n=n* / 


where 2 z? is the operator of the Heun equation such that (1) is written as JZH = 0. Besides, 

^ 2 ) ( z ) = - ^'( z ) + - (l—- + ip(z). 

J z z \ z z — 1 z — a J 

Obviously, ( 6 ) separates into two equations 


j?? [ Y C r, z 

\ n= 0 , 


Y s ^ n ) = 0 ’ 

\ 71 = 71 * / 

^ +Jf( Y SnA = 0 , 


( 6 ) 


(7) 

( 8 ) 


Now we collect in ( 6 ) terms having identical asymptotic nature as z —> 0. First we find that 
coefficients c n for n = 1,.n* — 1 are submitted to the recurrence (3): P n c n = Q n c n -\ + R n c n - 2 , 
where P n , Q n , R n are defined by (4) and the initial conditions are c_i = 0, cq = 1. 
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Further we find from (7) and ( 8 ) that the coefficients s n for n = n* + 1, n* + 2,... are submitted to 
the same recurrence relationships (3): P n s n = Q n s n -1 + R n s n - 2 , where s n »-i = 0 and another initial 
conditions include coefficients c n<t _i, c Uf - 2 : 

an* = c n ,-i [q- y(e + a5 - a - 1 )] - c nt _ 2 [(1 + 7 ) (2 - 8 - s) + a/3]. 

At the next step we can define coefficients c n for n = n* + 1, n* + 2,... From ( 8 ) we obtain the 
following relationship: 

PnCn — QnCn—1 A RnCn —2 “I - S n S n A T n S n —\ A U n S n —2, (9) 

where 


S n = a(l — 7 — 2 n), T n = £ + ad + (a + 1)(7 + 2n — 3), U n = 4 — 2n — a — (3. 

In this way for 7 £ {0} U ZA, using (5) we obtain a zero-exponent Frobenius solution, equal to 
unity at z = 0. The second solution locally can be defined as follows (see also (7)): 

OO 

Hs(a,q,a,/3,j,S-,z) = ^ s n z n , ( 10 ) 

n=n* 

where P n s n = Q n s n -i + Rn.s n -2 for n> n* and s nt = 1 , s nt _i = 0 . 

As it was mentioned above, c nt in (5) could be arbitrary. In other words, the linear combination 

HI (a, q, a, /3, 7 , 5; z) + C Hs(a , q, a, /3, 7 , 5; z ) 

is equal to unity at z = 0 for an arbitrary constant C. In this sense, for 7 £ {0} U ZA the choice of 
solution HI is non-unique. 

Consider now the function Hs(z) for arbitrary 7 . We should distinguish two situations: 7 = 1 and 
7 /I. For 7 / 1 we can use the following representation (see Table 2 in [9], index [0_][l+][a+][oo_]): 

Hs(a , q, a, (3, 7 , 6 ; z ) = z 1-7 HI (a, g - (7 — l)(e + ad),/3 - 7 + 1, a — 7 + 1, 2 - 7 , <5; z). (11) 

Notably, the later formula includes (10) as a particular case, justifying our way to introduce HI and 
Hs for non-positive integer 7 . 

For 7 = 1, repeating the arguments used to derive representation of HI in the case 7 £ {0} U Z~, 
we can find the following local representation 

OO OO 

Hs(a, q, a, f3 : 7 , <5; z) = ^ d n z n + log(z) ^ t n z n . (12) 

n =1 7i—0 

Here P n t n = Qntn-i + R n t n - 2 , where t -1 = 0 , to = 1 and (cf. (9)) 

Pndn = Qndn— 1 T Rndn—2 T S n t n T T n t n —1 A Urdn—2 • d—\ = do = 0 . 


4 Power series expansion at an arbitrary point 


Further we will extend the local Heun functions outside the circle of convergence of the series (2), (5), 
(12). To avoid numerical integration of the differential equation (1) we will use analytic continuation 
process and the power series expansion which is derived in this section. 

We seek the solution H( Zo t H' 0 )( a , q, a, P, 7, <5; z) to the equation (1) satisfying the conditions 


H(a, q, a, I3,i,d; z 0 ) = H 0 , 


d 

—H(a,q,a,fi,i,6;z ) 


Z=ZQ 


= < 


(13) 


Here zo is an arbitrary point that is assumed not to belong to the set {0,a, l,oo}. We will derive 
power series expansion of the Heun function in the following form: 


OO 

H (z 0 ,H 0 ,H' 0 ){a,q,Oi,P,i,8;z) = c, x (z - z 0 ) n . 

n=0 


(14) 
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Substituting (14) into (1) and collecting terms at the same power of z — zq, we obtain the following 
4-term recurrent relationship defining 47 


VnCn 


Qn.Cn —1 4“ T^n c n—2 4" S n Cn— 3; 


(15) 


where 


V n = -n(n - V)zq{z q - l)(z 0 - a), 

Qn = (n - 1) | [7 4 - 6 + e + 3(ra - 2 )] Zq + [(a 4 - 1 )(4 - 2 n - 7 ) - e — a<5] z 0 + 0(7 4 - n - 2 ) |, 

7 Z n = {(n - 2) [ 2(7 + 5 4- s) + 3(n - 3)] + otj3)zo - q - (n - 2) [(a 4- 1)(7 4- n - 3) + e + a5], 

S n = (n - 3) (7 4 - 6 + e + n - 4) + a/3. 

Obviously, the solution (14) satisfies the conditions (13) if the recurrence process starts with the initial 
conditions 

c-i = 0, co = Ho, ci = Hq. 

The series (14) converges inside the circle \z — zq\ < R, where R is the distance to the nearest 
singular point, R = min{|zo|, \zq — 1|, \zq — a|}. 


5 Basic algorithm 

First we consider evaluation of Hl(z) for 7 0 {0} U . We introduce the projection operator 
which, being applied to an analytic function, truncates its power series expansion at the point z to 
the first N terms. Using the expansion (2), we evaluate 

N N 

Hl)(z ) = ^ b n z n , (^Hl)'(z) = Y J ub n z n -\ (16) 

n =0 n=l 


as approximation of Hl(z ) and Hi!(z) in a vicinity of z = 0. 

In our algorithm the number N in the representations (16) is not fixed, it will be defined as we 
proceed with recurrent computation of b n (z) = b n z n and summation until a termination condition is 
satisfied. Namely, we stop the process when Hi) (z) and (^ 0 V_I Hlj (z) are not distinguishable 
in the used computer arithmetics and |&n(z)l < e > where e is the machine epsilon. 

To estimate the quality of the approximation, in view of (1) we compute the value 


Hl(z) = 


z(z— 1 )(z — a) (&o HI)”(z) + 7 (z — 1 )(z — a) + Sz(z — a) + ez(z — 1) Hl)\z) >, 


q — a/3z 

where HI)"(z) = Yln =2 n(n — 1) b n z n 2 . Then we suppose proximity of 


r 0 (z) = \Hl(z)-(^Hl)(z)\ (17) 

to the true error of the approximation | Hl(z) — Hl)(z)\. Note that numerical computation of 

Hl(z) near the point z = z* = q/(af3) is unreliable due to essential loss of significance. It can be 
suggested for a vicinity of z*, say, for {z : \q — ot/3z\ < 0 . 01 }, to use an estimate based on properties 
of the series, e.g. 

r 0 (z) = /N\b N (z)\ + eN\(^Hl)(z)\. (18) 

We can write the described algorithm as a function H[q (z) which returns 4-tuple 




where IV + 1 is the number of terms in power series, defined by the termination condition, r is the 
value computed with (17) or (18), / = Hl)(z) and f = Hi)'(z). 
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The scheme of computation of Hl(z) in the case 7 £ {0} U ZC is analogous, though a bit more 
involved. We use the expansion (5) and, instead of (16), define the function M'q(z) starting from the 
expression 

N N 

^2 c nZ n + log(z) ^2 S n zn - 

n=0, n^n* n=n* 

Assume that \z\ < >cRq, where Rq = min{l, |a|} is the radius of convergence of the series in (2) or 
(5) and x £ (0,1) is some coefficient chosen so that N defined by the termination condition can be 
expected to be moderate, say x = 0.5. Then we can use the numerical algorithm M/q for evaluation 
of the function Hl(z), its derivative and for estimation of the approximation error. 

Consider further the case \z\ > >cRq. First we define an auxiliary algorithm. Let zo be an arbitrary 
point not belonging to {0, l,a, 00}. Using (14) we define 

N N 

H (zo,Ho,%)){z) = ^Cn(z-Zo) n , (^2 H (z 0 , Ho ,H^y(z) = J2 n ~ • 

n —0 n— 1 

as approximations of H^ zoHqH ^{z) and Hq h')( z ) f° r z c l° se to z o • Here coefficients Cn are defined 
by the recurrence relationships (15). So we proceed with recurrent computation of c n (z) = Cn(z — zo) n 
and summation until the termination condition (analogous to that described above) is satisfied. 

We compute 


H 


(ZO ,H 0 ,H ’ 0 


)(z) = 


1 


q — a/3z 


z(z - 1 )(z - a)(&£ H ( Z o, Ho ,W))"(z) 


+ 


7 (z -l)(z-a) + Sz(z -a) + ez(z - 1) {^ 0 H {z 0 , Ho , H ' 0 ))'{z)\, 


and the value 

f(z 0 ,H 0 ,H^{z) = \H( ZOi H 0 ,H' q ){z) - H(z 0 , Ho ,H' q )) ( z )\, (19) 

supposing its proximity to the true error of the approximation \H^ zoHqH ^{z) — (^^ r .ff’( Z0i # 0i #'))(z)|. 
In view of essential loss of significance in computation of H^ Zq Hq h ^(z), near z = z* we define 

r( Z0 ,H 0 ,H^){z) = Vn\~c n (z)\ + eN\(^> z 2 H {z0jHo ^)(z)\. (20) 

We can write the described algorithm as a function !H^ zo Ho H ' q ){z) which returns 4-tuple 

% 0 ,H 0 ,H' 0 ) ■ z ^ [, f,f',r,N ], 

where N + 1 is the number of terms in power series defined by the termination condition, r is the 
value computed with (19) or (20), / = (^-H( Z0 ,.ff 0l tf'))(z) and f = (^>2 h (z 0 ,h 0 ,h^)\z)- 

Now we are in a position to describe the algorithm based on analytic continuation along a path 
from zero to z. Consider first the simplest version of the algorithm when the path is the line segment 
(0, z). At the first step, we compute 


[Hl 1 ,Hl' 1 ,r 1 ,N 1 \=M 0 (zi), 

where z\ = e iarg ^xi?o (see Fig. 2). Further, for p = 1, 2, and so on, we define 

R p = min{|z p |, | z p - 1|, \z p - a|}, 


Zp+l — 


z p + e iarg ^^xi? p if |z — z p \ > >cR p , 
z if | z — z p | < xR p , 


and compute 


[Hlp+li Hl p+i , T'p-i-i, IVp+l] — ?(( Zpj Hl p ,Hl'){Zp- |-l)- 
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Figure 2 : Analytic continuation using power series. 

The iterations stops when z p+ i = z. Finally, we have Hl p+ 1 , Hl' p+1 as approximations of Hl(z) and 
Hl'(z), respectively. We also compute the values = r\ +.. , + r p+ \ and Ny, = N± + ... + N p+ i+p+l. 
Here Ny, is the total number of power series terms which can be used as a measure of computer load 
and ry; may be an indicator of the approximation quality. 

The described above algorithm of continuation along a line segment is readily generalized for the 
case when 0 and 2 are connected by a polyline T. This gives us a way to compute the multi-valued 
Heun function. The resulting procedure can be considered as a function HI (T) which returns 4-tuple 

k-.T [/,/', r s ,IV s ], 

where / and f are the resulting approximations of the Heun function at z and its derivative. It should 
be noted that the function Hs(z ) for 7 / 1 is expressed through Hl(z) via (11) and to define Hs( T) 
through M'(Y) one should only compute the multi-valued function z -7 along T. In the case 7 = 1 the 
procedure of analytic continuation described above for HI can be applied with simple modification — it 
should start from the expansion ( 12 ). 

It is easy to observe that the size of the step in the described analytic continuation is small for 
points of the polyline T close to one of the singular points. This also means an increase of the number 
of used circular elements in the continuation procedure which, in its turn, may lead to loss of accuracy. 
The influence of the singular points can be reduced by a choice of the path of continuation. 

Let us give some details of application of the above algorithm for computation of the single-valued 
Heun functions. We will only use simple paths consisting of two or three line segments. We denote 
Ci = a and £2 = 1 if \a\ < 1 and £1 = 1 and C 2 = a otherwise. Let Rj = min{|£j|, |Ci ~ C 2 I} (. j = 1, 2 ) 
and rjj be the closest to Q point of the line segment ( 0 , 2 ) (see fig. 3). Note that Q / rjj due to the 
assumption z 0 &ioc U 3$ a oo- Let dj = \(j — rjj\ be the distance from (j to (0 ,z). If rjj is an internal 
point of (0, z) and dj < Rj/2 , then we introduce new point 

Sj = Cj + exp(i 7 r /2 + i arg (z)) min{Rj/2, \z - £j|} sign(Im(z/£j)). 

In this way we obtain the sequence of points n, which is either [0, z] or [0, <q, z\ or [0, z] or [0, < 7 , <,2 , z]. 
Then we define the path T that consequently connects the points of n by line segments (see an example 
of three-segment path in fig. 3) and define H£(z) = M'(T) and Hs(z) = T). 
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Figure 3: Path from zero to z consisting of three line segments 


6 Computation of single-valued Heun functions near singular points 

In this section we consider computation of single-valued Heun functions. We have already noted that 
the number of circular elements in the analytic continuation procedure becomes large as z approaches 
one of the singular points {l,o, oo}. So, we suggest an improvement of the algorithm; for this we 
will compute local solutions in a vicinity of the singular point, using results of [9] and writing these 
solutions in terms of Hl(z ) and Hs(z). 

Let z be close to 1. We use the following representation of the Heun function: 

HI (a. q, a , (3, 7 , 6; z) = C\ Hl(l — a, a/3 — q, a, f3, 5, 7 ; 1 — z) 4 - C 2 Hs( 1 — a, a/3 — q, a, /?, 5, 7 ; 1 — 2 ), 

( 21 ) 

where C \, C 2 are some constants, the Heun functions in the right-hand side (the first and the second 
local solutions in a vicinity of z = 1) correspond to the index [l + 0+][a+][oo+] in Table 2 [9]. There 
is a complication when the point 1 belongs to £$ a co (o E (0,1)). In this case the coefficients in 
( 21 ) are generally different for z belonging to the upper and the lower half-spaces: C^\ for 

{z : ± Im( 2 :) > 0 }. 

Consider branch cuts of Hl(l — a,...; 1 — z) and Hs( 1 — a ,...; 1 — z). The points of the branch cut, 
corresponding to SS aoo for HI (a ,...; z) and Hs(a ,...; z), are defined by the equation 1 — z = s(l — a), 
where s E (l,+oo). The set of points z = 1 + s(a — 1) constitutes the ray emanating from a to 
exp{iarg(a — l)}oo and, thus, it is located outside the circle of radius |a — 1| with center at z = l. 
Analogously, consider the branch cut of Hl(l — a ,...; 1 — z) and Hs(l — a ,...; 1 — z) corresponding 
to the branch cut BSioo- The points of the branch cut are defined by the equation 1 — z = s, where 
s E (1, + 00 ), thus, lying outside the circle of radius 1 with center at z = 1. The last possible branch 
cut of Hl( 1 — a ,...; 1 — z) and Hs{ 1 — a,...; 1 — z) is defined as 1 — 2 = q, q E (— 00 , 0) and coincides 
with the branch cut ^ioo of the function Hl(z ) in the left hand side of (21). 

In order to use the representation (21) we should define the coefficients C \, C 2 (or, in the special 
case, C| ±} , C^). Since for the general Heun equation an explicit solution to the two-point connection 
problem is not known, we will find the matching coefficients in the following way. First, we define a 
matching point z = z*, preferably being equally close to zero and 1 and distant from z = a. In our 
algorithm we fix 

*? = ! + -!! where »=/ 

2 v 2 I — sign(Im(a)) otherwise. 

Further, we apply the algorithms and Hs described in § 5 and find 

[fo,fo,r 0 ,N 0 ] = yff(a,g,a,/ 3 , 7 ,< 5 ,z*), 

[/i,/(,n, = MC(1 - a, a/3- q, a, /3,<5,7 , 1 - * 1 ), 

[/ 2 ,/ 2 ,r 2 , N 2 \ = Hs{ 1 -a, a/3- q,a,/3,5, 7,1 - z*). 
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Then we solve the linear system 


/o 

/o 



and obtain the matching coefficients C \, C 2 . It may also be reasonable to keep the found values 
C 1 = Ci (a, q , a, /?, 7 , 5), C 2 = C 2 (a, g, a, /3, 7 , (5) in computer memory. 

On finding Ci, C 2 (by computation or in the computer memory) we define 

:z^ [/,/',r,iV], 


where 


/ = C 1/1 + C 2 / 2 , f = -C 1 f 1 -C 2 f 2 , r = |Ci|n + |C 2 |r 2 , IV = 1^ + ^, 
[/i,/{,n,iVi] = #7(1 - a, a/3 - q,a,/3,8, 7 ,1 - 7 ), 

[/ 2 , f' 2 , r 2 , N 2 \ = Ms{ 1 - a, a/3 - g, a, /3, 8 , 7,1 - z). 

The described scheme can be repeated literally to define based on the representation 

Hs(a , g, a, /3, 7 , 5; z) = C[ Hl(\ — a, a(3 — q, a, (3,8 , 7 ; 1 — 2 ) + C 2 Hs(l — a, a/3 — g, a, /3, <5, 7 ; 1 — 2 ), (22) 

where C( and C 2 are some coefficients to be found (or two sets C’\ \ C' 2 ^ and C'\^\ C'^ —in the 
special case a £ ( 0 , 1 )). 

We note that finding Ci, C 2 or C(, C 2 demands computation of Hl(a ,...; z±) or Hs(a ,...; z±), 
and both terms Hl( 1 — a,...; 1 — z\) and Hs( 1 — a ,...; 1 — z±) in the right hand side of (21) or (22). 
So, if the matching constants are not known, the algorithms and are preferable over Jil 

and Ms in a sufficiently small vicinity of z = 1. In the code used in § 7 the algorithms are applied 
for \z — 1| < 0.05 i?i (R\ = min{l, |a — 11}) when the matching constants are not known and for 
\z — 1| < 0.25 R\ when the coefficients are already precomputed. 

In a very similar way we consider the case when z is located near a. We will use the following 
representation of the Heun function: 


Hl(a, q, a, /3, 7 , 8; z) 
= D 1 Hl 


a — 1 


, a/3 - -, a, /3, e, 7 ; 
a 


a — z 


+ D 2 Hs 


a -1 q a-z 

-, a/3 - a, p, e, 7 ;- 

a a a 


(23) 


where D \, D 2 are some constants and the two functions in the right-hand side correspond to index 
[a + 0 + l + ] [oo + ] in Table 2 [9]. There is a complication when the point a belongs to ^ioo or ^ooo- In 
this case we generally have different coefficients in ( 21 ) in the upper and the lower half-spaces: D^\ 
for {z : ± Irn(z) > 0}. 

Consider branch cuts of Hl((a — 1 )/a,... ;C) and Hs((a — 1 )/a ,where ( = (a — z)/a. The 
points of the branch cut that emanates from the singularity C, = (a — l)/a are defined in the z-plane 
by the equation a — z = -s(a — 1), where s £ (1, + 00 ). This set of points constitutes the ray emanating 
from 1 to exp{i arg(l — a)}oo and located outside the circle of the radius |1 — a| with center at z = a. 
The branch cut of Hl((a — 1)/a ,and Hs((a — 1 )/a,...; £) emanating from £ = 1 is a ray going 
from zero to exp{i arg(—a)}oo (z = a( 1 — s)); it is located outside the circle of radius |a| with center 
at z = a. The third possible branch cut of Hl((a — 1 )/a,...; £) and Hs((a — 1 )/a,...; £) is defined as 
a — z = aq , q £ (— 00 , 0) and coincides with the branch cut £$ a oo of the function Hl(a ,...; z). 

In order to find the coefficients D \, D 2 (or, in the special case, D^\ D^). first, we define a 
matching point z = z*, preferably being equally close to zero and a and distant from z = 1. We fix 


z 


* 

a 


a is 

— H-where s 

2 ^2 


sign(Im(z)) if a £ J? 0 oo U &ioo, 

a|a | -1 sign(Im(a)) otherwise. 
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Further, we find 


[fo,fo,ro,N 0 ] = Hi (a, q, a, [ 3 , 7, S, z*), 

[/ 2 ,/2,r 2 ,iV2] =*s 


a - 1 o Q a a ~ z a 

, G^p , G, (3 , £, p, 

a a a 

a - 1 0 <7 a - z* 

, Gp , g, p, p, 
a a a 


and D 1, Z) 2 are defined as solution to the linear system 


h 


h 


D 1 

-a“7i a-7P V^2 


/o\ 

fo)' 


It may also be reasonable to keep the found values D\ = D\(a, q, a, ( 3 , 7, 5 ), D2 = D-2(a,q,a, ,0,7, 5 ) 
in computer memory. 

On finding D\, D2 (by computation or in the computer memory) we define 


where 


^ (o) [/,/', r, N], 


f = Lfi/i + £>2/2, f = ~D 1 a- 1 f[-D 2 a- 1 f^ r = + \D 2 \r 2 , N = N 1 + N 2 , 


[f l ,f[,r 1 ,N l ] = XC 

lf2,f2,r2,N 2 }=rt 


a — 1 


a — z 


a /3 - 7, 

a a a 


a — 1 


a — z 


,a /3 - -,a,/ 3 , e, 7, 
a a a 


The described scheme can be repeated literally to define based on the representation 


Hs(a, q, a, / 3 , 7, < 5 ; z) 

= D\ HI 


a — 1 


a /3 -, a, ( 3 , e, 7: 

a 


a — z 


+ D'^Hs 


a — 1 


a — z 


,a /3 - -,a,/ 3 ,e,7; 

a a 


where and D 2 ( or ) dr the special case, D'^\ D' P) are some coefficients to be found. 

We note again that the algorithms Hf 1 ' 0 ' 1 and ‘Hs' 0 ' 1 are preferable over ‘M and Hs in a sufficiently 
small vicinity of z = a. In the code used in §7 the algorithms are applied for \z — a\ < 0.05 R a 
(R a = min{|a|, |1 — a|}) when the matching constants are not known and for \z — a| < 0.25 R a when 
the coefficients are already precomputed. 

Consider now vicinity of the point z = 00 . The branch cuts & 000 , &ioc and &aoo split the 
vicinity of infinity = {z : \z\ > max(l, |a|)} into three sectors, further denoted as S'd), S^ 2 \ S^ 3 \ 
Coefficients connecting the function Hl(z ) with two local solutions at infinity are defined for each of 
the sectors separately. Note that there is a special case when a is real, then the number of sectors 
decreases to two. 

We use the following representation of the Heun function near z = 00 for 2 : 6 S^ ( k = 1, 2,3): 


Hl(a , q , a, /3, 7 ,<5; z) = E^z a Hl ( ^ a ^ —@1 7 . a r £ _ ck,q: — ^ 7 + 1 , 0 : — 

\a a z / 

+ E^z~ a Hs(~, -—- —— + a(e — j3), a, a — 7 + 1, a — /3 + 1, d; -V (24) 
\a a z) 

(k) (k) 

where E\ , E\ are some constants. Further we will omit the superscript for brevity. The two 
functions in the right-hand side correspond to the index [oo + 0 + ][l + ][a_|_] in Table 2 [9]. 

Obviously, the branch cuts SS a oo and £%\oo transform for Hl(l/a,... ,1/z) and Hs(l/a,... ,1/z) 
to line segments connecting zero with z = a and z = 1 respectively. These branch cuts are located 
outside Coo. The branch cut &ooc suits both sides of (24). 

To find the coefficients E\, E 2 , first, we define a matching point z = z^. In our algorithm we fix 

zto = CooRooe iUk , for zeS^, k = 1,2,3, 
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where Uk is the angle of the mean line of the sector S^ k \ = max{l, |o|} and C 00 > 1 (= 2 in 
the computations, presented below). Then, we find 


[fo,fo,r 0 ,N 0 ] = HI(a,q,a,/3, 7,<5,z^), 

[fi, fi,n,N i] = HI f 1 , Q + —— + a(e - /3), a, a - 7 + 1, a - /3 + 1, S, , 

V a a z oo / 

[/ 2 , f 2 ,r 2 , N 2 \ = Hsf-, q + —— + a(s - /?), a, a - 7 + 1, a - /3 + 1, <5, , 

\ a a ^oo / 

and the matching coefficients are defined as solution to the linear system 

( \ pl\ = (fo\ 

\-[^w]“ a-1 (/l/«w + «/l) - [^] _a_1 (/2/Ac + “A)/ V-® 2 / VO/ 


Finally, we define a function returning values of f(z ) ~ Hl(z), f{z) « HI'{z) for sufficiently 

large |z|. On finding E \, (by computation or in the computer memory) we compute 




where 


/ = + E 2 z~ a f 2 , f = -E lZ - a -\f[/z + ah) - E 2 z~ a ~ 1 (f! 2 /z + af 2 ), 

r = \E lZ - a \n + \E 2 z~ a \r 2 , N = Ni + N 2 , 

[fi,fi,ri,Ni] = H((~, q + a ( S —@1 + a (e - f3),a,a - 7 + l,a - /3 + 1 ,( 5 , 

\a a z) 

[A, f 2 ,r 2 , N 2 \ = Hs ( 1 , q + —— + a(e - /3), a, a - 7 + 1, a - /3 + 1, 5, . 

\a a z) 

The scheme can be repeated literally to define Hs^°°^ based on the representation for z £ S^ 


Hs(a , q , a, /3, 7 ,5; z) = E'\ k ^z 01 Hl( -, At—i-— + a(e — j3),a , a — 7 + 1, a — (3 + 1, <5; -) 

V a a z / 

Q +—(_—A a / £ — /5 ),q;,ok — T + 1,ck — /3 + 1,<5; 

a z/ 


+ E’ (k) z~ a Hs[~, 


where E'^ and E'^ are some coefficients to be found. 

We note again that the algorithms yA°°) and Hs^°°^ are preferable over HI and Hs for sufficiently 
large \z\. In the code used in the following section the algorithms are applied for \z\ > 5 Coo Roo when 
the matching coefficients are not known and for \z\ > R oo otherwise. 


7 Numerical results 


In this section we present some results of numerical evaluation of the function HI. For this we will use 
the algorithm H[ with the improvements described in the previous section (HC^, Hl^ and HC^ 00 ^). 

It is known that in some cases Heun functions can be expressed in terms of hypergeometric functions 
(see [ 6 , 13, 8 , 16, 1] and references therein). For a test of the present algorithm we use the formula 
given in [ 1 ]: 


4 - A lb 2 ’*) -■ ( Hi 1 - <* - - 1) i - «*>» 


V4=z(l-z)’ 


where 2 Ei is the ordinary hypergeometric function. 

We define 

A( 2 ):=ffl( 4 ,^,!,1,2.7- 


y/4^(l -z) 
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Figure 4: Values of max{ —16,log 10 A(z)}. 
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and check numerically the identity A (z) = 0. Calculations are performed in the numerical computing 
environment GNU Octave and double precision (64-bit) arithmetics (the machine epsilon e is about 
2.22 ■ 1CT 16 ). 

Figure 4 shows in a semilogarithmic scale results of computations of the relative error 


|A(*)| | |A'(z 

i + \K z )\ i + W(z)\ 


on the grid (Re lm.z) G L(1000, [—20,20]) x L(1000, [—20,20]), where L(n,x ) is the set consisting 
of n linearly spaced in the interval x values (including interval’s end-points). It is easy to note that 
accuracy loss occurs mainly near the singularity z = 4. The maximum error for A(z) found on the grid 
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is 1.9635 • 1CT 14 . So, as one can see, the accuracy of computation of the function and its derivative is 
rather satisfactory for the used double-float arithmetics. 

In Fig. 5 we present (in a semilogarithmic scale) the value N(z) which means the total number 
of terms in power series used to compute Hl(A, |, |, g,2,z). This value varies from IV(0) = 1 and 
can characterize the time of computation. To be more specific, in the Table we present typical times 
of evaluation in Octave of Hl(A, |, |, |, 2, z) via the basic algorithm ‘HI and via the algorithm Hi 

improved by Hl^ l \ Hi^ and Hl^°°\ The computer in use has a 3.4 GHz Intel Core i5 CPU and 8 Gb 
of RAM. 


z 

Time (sec.), basic algorithm 

Time (sec.), 
1st evaluation 

improved algorithm 
subsequent evaluations 

20i 

0.04 

0.05 

0.007 

20 + ie 

0.075 

0.05 

0.007 

-20 

0.035 

0.05 

0.003 

0.99 

0.03 

0.02 

0.003 

4 + 0.01i 

0.075 

0.05 

0.002 
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